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Abstract 

We address via numerical simulation the two-dimensional bimolecular annihilation 
reaction A+A — > in the presence of quenched, random impurities. Renormalization 
group calculations have suggested that this reaction displays anomalous kinetics at 
long times, c&{t) ~ at 6 ^ 1 , for certain types of topological or charged species and 
impurities. Both the exponent and the prefactor depend on the strength of disorder. 
The decay exponents determined from our simulations agree well with the values 
predicted by theory. The observed renormalization of the prefactor also agrees well 
with the values predicted by theory. 

1 Introduction 



Surface reactions in two-dimensions show a variety of interesting behavior. 
Qualitatively, the difference between two and three dimensions is that diffu- 
sive mixing is less effective in two (and one) dimensions. With less mixing, 
transport limitations become more significant, and this leads to a breakdown 
of the law of mass action in two dimensions. Formally, two-dimensions is the 
upper critical dimension for bimolecular surface reactions [1,2]. The law of 
mass action, or local chemical kinetics, is the simplest theory with which to 
predict the reactant concentration in the A + A —> reaction: 



IT ~ ~ k ° A 

c A (0) = n , (1) 

where k is the reaction rate constant. This theory fails because, without input 
reactants, the reaction is diffusion limited at long times. The approach of 
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Eq. (1) does not take into account diffusion limitations, and so this theory 
cannot be a predictive one. If, in addition, the reaction occurs in a disordered 
medium, there will be an interplay between the effects of disorder and diffusion 
limitations, and the simple law of mass action can make no statements about 
the resulting kinetics. Finally, the law of mass action cannot account for any 
clustering or other collective behavior of the reactants, since it effectively 
assumes that the reactants are perfectly well mixed. 

A more sophisticated level of theory for chemical reactions is the reaction 
diffusion partial differential equation: 



(If \ 

= DV 2 c A + (3DV ■ (c A Vv) - kc\ . (2) 

We will be interested in the case of random initial conditions, where ca(x, 0) 
is a Poisson random number with average no- Here D is the diffusivity, v(r) is 
the quenched potential in which the reactants diffuse, and (3 = l/{k-Q,T) is the 
inverse temperature. This approach, surprisingly, is also qualitatively wrong 
in two dimensions. This approach fails because it is a type of mean field 
theory, and two dimensions is the upper critical dimension for bimolecular 
kinetics, below which mean field theory fails. If one insists to use Eq. (2), 
one could predict the true behavior only by using a renormalized, effective 
reaction rate in place of the bare rate. This effective reaction rate is exactly 
what the renormalization group treatment of this problem predicts [1,3-5]. 
What does the reaction diffusion equation leave out, which causes it to miss 
the renormalization of the effective reaction rate? Technically this approach 
does not capture the integral occupation number constraint at each site. That 
is, at each place in space, there can be only an integer number of reactants. 
This seemingly trivial constraint is what leads to a renormalization of the 
effective reaction rate. 

Previous simulation work has been conducted for a variety of simple reactions 
in two dimensions. Many simulations have been performed on clean metal 
surfaces (for a review see [6]). Simple systems, such as oxidation of CO on 
single crystal Pt(110), show surprisingly rich behavior [7], ranging from spi- 
rals and standing waves to chemical turbulence [8,9]. Simulation studies have 
confirmed several analytical results for the simpler systems (see, for example, 
[10]). Fascinating effects have been observed on surfaces with patterned dis- 
order [11,12]. Surfaces with spatially random adsorption energies have been 
shown to lead to a variety of phases observed at steady state [13]. Finally, 
simulations have been extended to the case of fractal media (see, for example, 
[14,15]). 

Unlike these previous simulation studies, we are interested in the kinetics of 
the ionic reaction A + + EP — > in two dimensions in the presence of quenched, 
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charged disorder. The disorder, and the reactants themselves, can either be 
simple electrostatic charges or topological "charges" such as dislocations or 
disclinations. For either system, the quenched defects interact with the diffus- 
ing reactants via a long-ranged, Coulombic potential energy. The upper critical 
dimension for this diffusive process is also two. The effects of this disorder are 
dramatic, and, in fact, the potential energy landscape induced by the random 
charges is so rugged that the motion is sub-diffusive in two dimensions [16,17]. 
In other words, this special type of disorder causes the mean square displace- 
ment of random walkers to increase sub-linearly with time, (r 2 (t)) ~ bt 1 ^ 5 . 
The scaling exponent is continuously variable in the strength of disorder and 
can be found exactly [18-24]. 

Each of the reactions A + A ^ 0, A + B ^ 0, and A + + B # AB have 
been studied in the presence of quenched, singular disorder by field-theoretic 
methods [4,5,25]. Renormalization group theory was used to derive the long- 
time values for the reactant concentrations. Interestingly, the A + + B # AB 
ionic reaction at high temperature behaves like the A + A — > model when 
the back reaction rate is equal to zero [25,26]. This is because the Coulomb 
interaction between the reactants inhibits the reactant segregation that would 
otherwise occur in this reaction. We chose to study the A + A — > reaction 
as a simplified model. 

We can understand the results of these renormalization group studies by a 
simple argument. The A + A — > reaction becomes diffusion limited at 
long times in the quenched disorder. Simple theory for diffusion-limited re- 
actions, therefore, predicts that the effective rate constant will be propor- 
tional to the diffusivity. Since the diffusion is anomalous in this disorder, 
(r 2 (t)) ~ ADt(t/t y s , we should expect k cS oc D cS ~ D(t/t y s , where 5 
is the same measure of the strength of disorder as above. Here D is the bare 
diffusivity, and t = (Ar) 2 /(2D) is a characteristic time for diffusion on a lat- 
tice of spacing Ar. We then might expect to use mean field theory with this 
effective reaction rate to predict the concentration of reactants 



where k* is some fixed-point effective reaction rate to be determined by a 
more careful calculation. This simple argument is in agreement with rigorous 
renormalization group results, which predict [4] 



We here test these renormalization group predictions. We analyze via numer- 
ical simulation the behavior of the A + A — > reaction in two- dimensions in 




as t 



oo 



(3) 



k* = 3p 2 -fD . 



(4) 
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the presence of singular disorder. We calculate by simulation the concentra- 
tion profiles at long times for various strengths of disorder. We compare both 
the exponent and the prefactor predicted by renormalization group arguments 
against those measured via simulation [4]. This paper is organized as follows. 
In Section 2 we describe the disorder that is appropriate for ionic systems. We 
express the effects of the quenched disorder in terms of a random potential. 
In Section 3, we describe the master equation that defines our model of the 
reaction. We describe a method for numerically solving the master equation 
by a Poisson process that is effective for strong disorder. In Section 4, we take 
an alternative approach and describe an exact stochastic partial differential 
equation (SPDE) that can be derived by field-theoretic means from the mas- 
ter equation. We describe a method for solving the SPDE numerically that 
is effective for weak disorder. In Section 5, we implement these solution tech- 
niques and analyze the long-time behavior of the reactant concentration at 
various strengths of disorder. We discussion our numerical results in light of 
the theoretical predictions in Section 6. We conclude in Section 7. 



2 The Disorder 

The two-dimensional ionic reaction A + + B — > behaves like the nonionic 
bimolecular annihilation reaction A + A — > because the Coulomb interaction 
inhibits the reactant segregation that naturally occurs in the A + B — > 
reaction [25,26]. This, then, implies that the reactant concentration for the 
A + + B ^ model is proportional to that for the A + A — > model at 
long times. The presence of additional, quenched, charged disorder does not 
change this argument. The A + A — > reaction is much easier to simulate than 
the ionic A + + B — > reaction, since one does not have to track the long- 
range forces between the reactants. For this reason, we choose to simulate the 
A + A — > model as a general test of the renormalization group predictions 
for two-dimensional bimolecular reactions in the presence of singular disorder 
[4,5,25]. 

The type of physical systems that we have in mind are systems with "ionic" 
disorders. Such disorder could occur, for example, on the surface of an ionic 
crystalline lattice. Imagine that the lattice has line defects caused by dislo- 
cation line pairs forming line vacancies or line interstitials. If these defects 
were immobile, they would generate a random, quenched electrostatic poten- 
tial on the surface. The interaction between ionic reactants on the surface and 
these line charges would be logarithmic, which is technically "relevant" in two 
dimensions. The 1/r interactions between the ions themselves is technically 
"irrelevant" and can be ignored. 

Alternatively, our system can be viewed as modeling the dynamics of a collec- 
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tion of line dislocations in a solid. The hexatic-solid transition in two dimen- 
sional monolayers occurs by pairing of dislocations and is a good example of an 
ionic reaction, albeit with vector "charges" [27]. If there were some additional 
dislocations pinned by defects, these defects would interact logarithmically 
with the moving dislocations, just as in our model. A physical system that 
more closely resembles our model is the liquid-hexatic transition in two di- 
mensional fluids. Here the transition occurs by pairing of disclinations, which 
are a perfect manifestation of topological charges that interact logarithmically 
[27]. Additional, pinned disclinations would generate the quenched disorder 
present in our model. 

Experiments have, indeed, seen effects of pinned disorder on the liquid-hexatic 
transition for both flux lines in charged density wave systems [28,29] and 
surfactants in hexatic multilayers in Langmuir-Blodgett films [30]. Detailed 
renormalization group calculations are in agreement with the observed effects 
[25] . An additional physical system, for which experiments could be performed, 
are charged vortices in a thin fluid film between two spatially-addressable 
electrodes [31]. 

We use a quenched random potential to represent the potential induced by 
the immobile defects. The random potential is Gaussian at long wavelengths 
[32] , and the appropriate form of the potential-potential correlation function 
at long wavelengths is 

* W (*0~7A 2 - (5) 



Here 7 is a measure of the strength of disorder. It is roughly proportional to 
the square of the density of defects. Our simulation takes place on a square 
lattice, and a natural form of the potential-potential correlation function is 

* vv{kx > ky) = 4-2cos(fc x Ar)-2cos(jfc I/ Ar) ' (6) 



where Ar is the lattice spacing. This form reduces to that of Eq. (5) at long 
wavelengths, k — > 0, and in the limit of small lattice spacing, Ar — > 0. We 
generate the random potential v(r) by first generating -O(k) in Fourier space, 
where the values at different wave vectors are independent Gaussian random 
numbers, and then performing an inverse fast Fourier transform [33,34]. We 
use periodic boundary conditions for the random lattice. 

A random potential of this form is known to have a significant effect on dif- 
fusing species, resulting in an anomalous mean square displacement at long 
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times [16-24] 



(7) 



The scaling exponent, 1 — 5, depends on the strength of disorder. For Gaussian 
disorder, the exponent is exactly given by 



5 = 



1 + 



57T 



P 2 1 



(8) 



3 The A + A -> Reaction 



Our model of the A + A — > reaction can be defined mathematically by a 
master equation. We assume that the reaction takes place on a square N x N 
lattice with lattice spacing Ar. The reactants move diffusively in a random 
potential on the lattice. The reactants annihilate at rate k whenever two meet 
on the same lattice site. The master equation describes how the probability 
of each possible configuration of particles on the lattice, P({rii},t), changes 
with time. The appropriate master equation for our system is 



dP({n,},t) =E[r _ lK + l)p{ ni _ ljn . + x _ T -, niP] 

ij 

+ 2(ihf + 2)(Ui + 1)p( -' Hi + 2 ' t] 

-n t (n t -l)P] (9) 

where rii is the number of particles at site i, D is the diffusivity, k is the reaction 
rate constant, and Ar is the lattice spacing. Here r^" 1 , to be defined below, is 
the rate at which particles hop from lattice site % to nearest neighbor lattice 
site j. The summation over i is over all sites on the lattice and the summation 
over j is over all nearest neighbors of site i. We place the reactants randomly 
on the lattice at time t = 0. The initial concentration at any given site will, 
thus, be Poisson, with average density that we define to be no: 



nil 

(n t /(Arf)=n , (10) 



This master equation can be solved exactly by directly considering the particle 
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reaction and diffusion process [35]. The master equation implies that a given 
configuration of reactants on the lattice follows a continuous-time, Markovian, 
Poisson process. In this process, all possible changes to the lattice configuration 
occur at specific rates. That is, both the reaction and diffusion moves have 
rates. The rates of diffusion depend on the values of the potential on the 
lattice. Certain moves are favored over others, due to local gradients in the 
potential. The rate of hopping to a nearest neighbor site is r^ 1 . There are four 
of these rates at each lattice site. The rate for hopping to the right is 

r R\ x ,y) = j^^ e w{P[v(x,y) ~v(x + Ar,y)]/2} , (11) 



for hopping to the left is 

Tl\x, y) = exp {p [v(x, y) - v(x - Ar, y)] /2} , (12) 



for hopping above is 



T v \x,y) = J^y^viPHx^y) -v(x,y + Ar)]/2} , (13) 



and for hopping below is 



t d 1 V) = J^yl ex P if 3 H x i y)~v(x,y- Ar)] /2} . (14) 



At each lattice site, there are rii{rii — l)/2 possible reaction moves. Each of 
these moves has the rate 



r rxn (Ar) 2 ' 



Overall, for the whole lattice there are 4J]i^i possible diffusion moves and 
J2i n i{ n i — l)/2 possible reaction moves. 

A Markov process is initiated by placing the particles randomly on the surface, 
with average density hq. The Markov chain is then generated by picking one 
of these possible events according to its probability of occurrence and incre- 
menting time according to its rate of occurrence. The probability of event a 
occurring, out of all the possible diffusion and reaction moves, is 

r- 1 

P(event a) = a — . (16) 



7 



Since the process is Poisson, time is incremented by 



where ( is a uniformly distributed random number between zero and one. This 
step of the Markov chain is repeated until only zero or one particles remain 
on the lattice. 

We find this procedure to be efficient when the strength of disorder, /3 2r y, 
is large. For computational convenience, we make two approximations. First, 
when a particle moves to a site that is already occupied by another particle, the 
annihilation reaction occurs immediately. In other words, we set the reaction 
rate, k, equal to infinity. Note that the renormalization group predictions 
are independent of the bare reaction rate [4]. Physically this is because the 
reaction is diffusion-limited at long times, and so the bare reaction rate does 
not matter. Second, we assume that each particle has an equal probability of 
being moved 

P(moving particle a) — — , (18) 

n 



where n is the total number of particles on the lattice. The direction of move- 
ment of a chosen particle, however, is specified by the local gradient of the ran- 
dom potential, as in Eqs. (11-14). Specifically, we generate another uniformly 
distributed random number, which is compared against the four different hop 
probabilities: 

P(hop i) = • (19) 



The chosen move is performed, and time is incremented by At = N 2 / {n T ij 1 ) 
The uniform choice of reactant to move and the approximate time incremen- 
tation are both exact in the long-time limit when Ar — > in the hopping rates 
(11)-(14). 



4 The Stochastic Equation 



Alternately, we can derive a stochastic partial differential equation by mapping 
the master equation onto a field theory. The field theory looks like a Bosonic 
quantum field theory due to the constraint of an integral occupation number at 
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each lattice site. The continuous-time stochastic partial differential equation 
for the A + A — > reaction is [4] 



del 

= DV 2 a + (3DV ■ (aVu) - ka 2 + irja 
a(x, 0) =n , 



(20) 



where the real, Gaussian random field r}(x.,t) has zero mean and variance 
(?y(x, t)i](x.', t')) = k5(~x — x.')5(t — t'). The physical concentration is given by 
taking the average of the solution over the random field r\. Since the reaction 
occurs on a lattice, the spatial derivative symbols are actually a shorthand 
for finite difference expressions. The field-theoretic derivation shows that this 
stochastic differential equation should be interpreted in the ltd sense. 

The stochastic partial differential equation was solved numerically by using a 
second-order, semi-implicit method. This method gives exact solutions to the 
master equation without any approximations, in the limit that the integration 
time step tends to zero. The process of solving the stochastic partial differential 
equation takes place on a square lattice. We define b[a(x,y,t)} = D'V 2 a — ka 2 
and [i = y/k/ Ar. An explicit method for Eq. (20) equation is given by [36] 



a(x, y,t + At) = a(x, y, t) 
+ 



^(b[a(x,y,t)] +b 



2 I 



a(x, y, t) + i(iy/ia(x, y, t)z(x, y, t) 



+ Atb[a(x,y,t)} 



At 



a(x,y,t) + ^b[a(x,y,t)} 



z(x,y,t) 



H a(x,y,t)[z(x,y,t) - 1] 



(21) 



Here z(x,y,t) is a Gaussian random number of unit variance. These random 
numbers are uncorrelated at different positions and times. This numerical 
method is accurate to 0[(At) 2 } in the weak sense. A semi-implicit method is 



on the right hand side with 
so accurate to 0[(At) 2 ] in the 



generated by replacing b a + \\i\ftaz + Atb{a 
b[a(x,y,t + At)}. This semi-implicit method is a 
weak sense. Direct implementation of this approach leads to a complicated, 
non-linear matrix equation to solve. We use, instead, an operator-splitting 
approach [37]. On even integration time steps, we use a method implicit in 
the x-diffusion terms and explicit in the y-diffusion and reaction terms. On 
odd integration time steps, we use a method implicit in the y-diffusion terms 
and explicit in the x-diffusion and reaction terms. This operator-splitting, 
semi-implicit approach leads to tridiagonal matrix equations that are simple 
to solve. We found that including a time step implicit in the reaction terms was 
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Fig. 1. The reactant concentration as a function of time for (3 2 ^ = 1, k = 1, 
and no = 1. The simulation data are shown by the solid line. The long-time value 
predicted by renormalization group arguments is shown by the long-dashed line. 
The concentration predicted by renormalization group calculations for all times is 
shown by the short-dashed line. At times longer than those shown here, the two 
dashed lines converge. 

not useful. The periodic boundary conditions enter these equations through 
the finite-difference definition of the diffusion terms. 



5 Results 



For weak disorder, (3 2 ^ < 1, the master equation was solved by numerically 
integrating the exact stochastic differential equation. The concentrations pro- 
duced in a typical simulation run are illustrated in Figure 1. The simulation 
values at short times are quite far off from the long-time values predicted 
by the renormalization group treatment. This is because the asymptotic scal- 
ing of the reaction concentration is reached only for very long times for the 
parameters of Figure 1. At very long times, for a suficientiy large lattice, the 
concentration observed by simulation would be the asymptotic value predicted 
by theory. To achieve an agreement between theory and simulation, we also 
show in Figure 1 the concentration predicted for all times with the use of the 
running coupling in Eq. (1). This approximate concentration profile is derived 
in a straight-forward fashion from the flow equations for this system [4]: 



c{t) = W^tl 

U n (t/t o y-i + k(l*)t 
k(n 3 ^ D (22) 

As we can see from this more detailed theoretical prediction, the renormal- 
ization of the effective reaction rate to the asymptotic value is rather slow, 
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occurring as (t/t )~ 3S . At long times, the level of noise in the simulation data 
becomes more significant. This is because there are fewer particles at longer 
times, and so the statistical error is greater. If a greater lattice size or a smaller 
time step were used, then less noise would be observed at these times. Bal- 
ancing these considerations against computational feasibility, we chose to use 
2048 x 2048 square lattices. 

The slow renormalization of the effective reaction rate prevents a direct sim- 
ulation of the reaction for arbitrary values of the parameters. To counter- 
act the effect of the slow renormalization, we performed additional simula- 
tions by starting with the renormalized reaction rate predicted from theory, 
k = k* = 3/3 2 7-D. For runs with this initial value of the reaction rate, the 
concentrations reach the asymptotic scaling at fairly short times. 

To determine the observed values of the the decay exponent, 1 — 5*, and the 
renormalized reaction rate, k*, a log-log plot is made for the concentration of 
species A versus time. The slope of this plot gives us the renormalized decay 
exponent, and the y-intercept gives us the renormalized reaction rate. With 
this method, far more data points will be plotted and used at long times than 
intermediate times in the fit to determine the renormalized values. Therefore, 
the values determined from the fit will be heavily biased by data points gath- 
ered at long times. This biasing is artificial, since the concentration profile 
must be continuous as a function of time, and so we should not overweight 
the long-time section. To eliminate this bias, we used data points in powers 
of two to perform the fits. For example, data points corresponding to times 
steps 1, 2, 4, 8, 16, 32, 64, . . . were used to determine the slope and the pref- 
actor. Note that to determine k* we need to assume a value for to- We use 
to = (Ar) 2 /(2£>), which is suggested by simple perturbation theory for this 
problem. 

Figure 2 compares the exponent of the concentration decay observed by simu- 
lation with that predicted by theory. Figure 3 compares the effective reaction 
rate observed by simulation with that predicted by theory. These two plots 
show that the theoretical values agree well with those obtained by simulation. 
The error bars were determined by considering both random and systematic 
errors. We found that the most conservative method was to vary the smallest 
time at which we start the fit. We found that the slope and prefactor were 
somewhat sensitive to which data points were used in the calculations. The 
error bars encompass both the random, statistical errors and the systematic 
error associated with the range of data that we fit. 

For strong disorder, f3 2r j > 1, the approximate Poisson process was used to 
solve the master equation. For each set of parameters, we performed 10 simu- 
lations on 2048 x 2048 lattices to generate average concentration profiles. As 
before, a logCA versus logt plot of the simulation data was used to determine 
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Fig. 2. Presented are the observed values for the exponent of the reaction con- 
centration decay, determined from ca(£) ~ at 6 *" 1 , as a function of the strength of 
disorder. Renormalization group predictions are shown by the solid line. 
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Fig. 3. Presented are the observed values of the renormalized reaction rate, de- 
termined from ca(£) ~ (t/to) s * /(k*t), as a function of the strength of disorder. 
Renormalization group predictions are shown by the solid line. 

the decay exponent, 1 — 5*, and the renormalized reaction rate, k*. Also as 
before, we used data exponentially spaced in time. For strong disorder, the 
renormalization of the effective reaction rate constant occurs very quickly. 

Figure 4 compares the exponent of the concentration decay observed by simu- 
lation with that predicted by theory. Figure 5 compares the effective reaction 
rate observed by simulation with that predicted by theory. 



6 Discussion 

The most significant prediction of the renormalization group studies is the 
decay exponent, c\(t) ~ at 1 ~ s , where 5 is given by Eq. (8). All simulations 
for weak disorder, f3 2r y < 1, show an observed slope that is in agreement 
with this prediction, to within the error bars. This is significant, since the 
renormalization group studies are, in principle, expansions in the parameter 
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Fig. 4. Presented are the observed values for the exponent of the reaction con- 
centration decay, determined from ca(£) ~ at 6 *" 1 , as a function of the strength of 
disorder. Renormalization group predictions are shown by the solid line. 
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Fig. 5. Presented are the observed values of the renormalized reaction rate, de- 
termined from ca(£) ~ (t/to) s * /(k*t), as a function of the strength of disorder. 
Renormalization group predictions are shown by the solid line. 

/9 2 7, and so the predictions are strictly accurate only for weak disorder. It was 
anticipated in the renormalization group treatment, however, that the result 
for the decay exponent may be exact to all orders of (3 2/ y [4] . 

For strong disorder, f3 2r j > 1, some small discrepancies exist between the simu- 
lation and the theoretical predictions for the decay exponent. The approximate 
Poisson process method that we used to solve the master equation requires 
the effective reaction rate to renormalize from infinity to k*. This process oc- 
curs slowly for weak disorder, as shown by Eq. (22), and this is the reason for 
the discrepancy for 2r y pa 1. Simulations on lattices larger than 2048 x 2048, 
which would provide data to lower concentrations, would eventually yield an 
exponent in agreement with the predicted value. 

The renormalization group studies also predict the fixed-point value of the 
effective reaction rate, Eq. (4). This prediction, strictly speaking, is an ex- 
pansion in the disorder strength, and we would generally expect there to be 
corrections higher order in j3 2r y. All of the simulations for weak disorder are 
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in agreement with this prediction. For (3 2/ ~f m 0.33, there is a small discrep- 
ancy not contained within the error bars. The numerical integration of Eq. 
(20) is exceedingly difficult for large j3 2r y, and the small discrepancy between 
simulation and theory is most likely due to use of a time step that was not 
quite small enough. Resolution of this discrepancy would have required an 
integration time step substantially smaller than is feasible, even on our high- 
performance workstations. 

Even for strong disorder, the simulations show an effective reaction rate in 
agreement with the predicted values. This is significant on two counts. First, 
the renormalization group predictions are only a first order approximation 
for k*, and it is surprising that the higher order corrections are so small for 
/3 2 7 ~ 20. Second, the simulation technique for strong disorder requires the 
effective reaction rate to renormalize from infinity to a finite k*. That this 
happens so effectively is impressive. The discrepancies in the observed effective 
reaction rate for f3 2r y ~ 1 are, again, due to the slow renormalization of the 
effective reaction rate. 



7 Conclusions 



Our numerical simulation of the A + A — ► reaction is in agreement with 

the field-theoretic renormalization group predictions [4]. An anomalous de- 
cay exponent, predicted by the renormalization group studies, is observed. No 
significant discrepancy between the numerical and analytical predictions are 
observed over a wide range of disorder strengths. More impressively, the effec- 
tive reaction rate observed from the numerical simulations is in quantitative 
agreement with the renormalization group predictions for the same wide range 
of disorder strengths. This agreement is unexpected and may signify that the 
one-loop renormalization group results are more accurate than anticipated. 
With these simulations, we now have satisfying agreement between rigorous 
field-theoretic results, simple physical arguments, and exact numerical results 
for this interesting case of anomalous kinetics. 
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